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We give a rigorous and mathematically clear presentation of the Covariant and Gauge Invariant 
theory of gravitational waves in a perturbed Friedmann-Lemaitre-Robertson- Walker universe for 
Fourth Order Gravity, where the matter is described by a perfect fluid with a barotropic equation 
of state. As an example of a consistent analysis of tensor perturbations in Fourth Order Gravity, we 
apply the formalism to a simple background solution of R" gravity. We obtain the exact solutions 
of the perturbation equations for scales much bigger than and smaller than the Hubble radius. It 
I is shown that the evolution of tensor modes is highly sensitive to the choice of n and an interesting 

, new feature arises. During the radiation dominated era, their exist a growing tensor perturbation 

for nearly all choices of n. This occurs even when the background model is undergoing accelerated 
^ ' expansion as opposed to the case of General Relativity. Consequently, cosmological gravitational 

wave modes can in principle provide a strong constraint on the theory of gravity independent of 
other cosmological data sets. 
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I. INTRODUCTION 



"Y* '' In the near future. Gravitational Waves (GW) will become a very important source of data in cosmology. Cosmo- 
5-H logical GW are produced at very early times in the evolution of the universe and almost immediately decouple from 
the cosmic fluid. Consequently they carry information about the conditions that existed at this time, thus providing 
a way of constraining models of inflation [Ij . 
CNJ Even if GW are decoupled from the cosmic fluid their presence still influences some features of the observable 

^ universe. In particular, a GW background will produce a signature that can be found in the anisotropics ^ and 
polarization of the Cosmic Microwaves Background (CMB). This, together with the remarkable improvements 
^0 in the sensitivity of CMB measurements, opens the possibility of obtaining important information about GW in an 
indirect way. 

In the last few years the idea of a geometrical origin for Dark Energy (DE) i.e. the connection between DE and a 
non-standard behavior of gravitation on cosmological scales has attracted a considerable amount of interest. 

Higher order gravity, and in particular fourth order gravity, has been widely studied in the case of the Friedmann- 
Lemaitre-Robertson- Walker (FLRW) metric using a number of different techniques (see for example [1, [R [E, 0J1J1, 
^d\). Recently a general approach was developed to analyze the phase space of the fourth order cosmologies[lll [l3l[l4| . 
providing for the first time a way of obtaining exact solutions together with their stability and a general idea of the 
qualitative behavior of these cosmological models. 
\^ ' The phase space analysis shows that for FLRW models there exist classes of fourth order theories in which the 
5^ I cosmology evolves naturally towards an accelerated expansion phase which can be associated with a DE-like era. 
Although this feature is particularly attractive, a problem connected with the use of these theories is that there is too 
much freedom in the form of the theory itself. Consequently, it is crucial to investigate these models in some detail 
in order to devise observational constraints which are able to eliminate this degeneracy. 

A key step in this process is the development of a full theory of cosmological perturbations. A detailed analysis 
of the evolution of the scalar perturbations on large scales has recently been given in [ISj . Here we will focus on the 
evolution of the tensor perturbations, which are related to GWs. This is motivated by the well known fact [l^, that 
the features of GWs in GR are rather special and therefore the detection of any deviation from this behavior would 
be a genuine proof of the break down of standard GR. 

The aim of this paper is to present a general framework within which to consistently analyze tensor perturbations 
of FLRW models in fourth order gravity (see [i1[i1[i1[1q| for other recent contributions to this area). As an explicit 
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example we apply our approach to the case of i?" gravity. We investigate the possible constraints one can place on 
such a model through future observations of gravitational waves independently of existing cosmological data. 

In order to achieve this goal a perturbation formalism needs to be chosen that is best suited for this task. One 
possible choice is the Bardeen metric based approach f2l|, [S^, |23'| which guarantees the gauge invariance of the results. 
However this approach has the drawback of introducing variables which only have a clear physical meaning in certain 
gauges [13]. Although this is not a big problem in the context of General Relativity (GR), this is not necessarily true 
in the case of higher order gravity and consequently can lead to a miss-interpretation of the results. 

In what follows we will use, instead, the covariant and gauge invariant approach developed for GR in [2^. [25l. [26l. [27l. 
[28l . [29j which has the advantage of using perturbation variables with a clear geometrical and physical interpretation. 
We take advantage of the fact that in this approach the non-Einstein part of the gravitational interaction can be 
considered as an effective fluid (the curvature fluid) coupled with standard matter. This specific recasting of the field 
equations makes the development of cosmological perturbation theory even more transparent. 

The main results of the paper are as follows. (1) We find that the evolution of tensor modes is extremely sensitive 
to the choice of f{R) theory. (2) In the specific case of i?" gravity, the tensor modes are in general weaker due to a 
higher expansion rate in the background. (3) During the radiation dominated era, their exist a growing tensor mode 
for nearly all interesting values of n. 

The paper is organized as follows. In section II we give a brief review of the H-3 gauge invariant covariant approach 
in a general setting. In section III we present the equations necessary for the study of linear tensor perturbations for a 
general imperfect fluid. In section IV we investigate how these equations are modified when considering fourth order 
gravity. In section V we adapt these equations for the specific case of i?" gravity and study tensor perturbations both 
in vacuum and in the presences of dust/radiation fiuid. Finally, we present our discussions and conclusions in section 



The starting point (and the corner stone) of our analysis is the 1-1-3 covariant approach to cosmology ^3Q] . This 
approach consists of deriving a set of first order differential equations and constraints for some suitable, geometrically 
well defined quantities (the 1+3 equations) that are completely equivalent to the Einstein field equations. This has 
the advantage of simplifying the analysis of general spacetimes which can be foliated as a set of three dimensional 
(spacelike) surfaces. In the following we give a very brief introduction to the parts of this formalism used in this 
paper. 



We will adopt natural units {h — c = = SttG = 1) throughout this paper, Latin indices run from to 3. The 
symbol V represents the usual covariant derivative and d corresponds to partial differentiation. We use the — , -f , -I- 
signature and the Riemann tensor is defined by 



VI. 



II. THE 1+3 COVARIANT APPROACH TO COSMOLOGY 



A. 



Preliminaries 



where the W^bd is the usual Christoffel symbol (i.e. symmetric in the lower indices), defined by 

W'bd = y'"' {gteM + 9edM - 9bd,e) ■ 



(1) 



(2) 



The Ricci tensor is obtained by contracting the first and the third indices 



Rab 



g'^ Rcadb ■ 



(3) 



Finally the Einstein-Hilbert action in the presence of matter is defined as 




(4) 



B. 



Kinematics 



In order to derive the H-3 equations we have to choose a set of observers i.e. a 4-velocity field u". This choice 
depends strictly on the theory of gravity that we are treating. In this section we give the set of equations for a general 
velocity field. In later sections we will discuss how this situation is modified in the case of f(R) gravity. 
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Given the velocity u"", we can define the projection tensor into the tangent 3-spaces orthogonal to the flow vector: 

Kb = gab + UaUb => h°-hh^c = /l"c , habu'' = 0, (5) 

and the kinematical quantities can be obtained by splitting the covariant derivative of Ua into its irreducible parts: 

VfaUa = WbUa - AaUb , V bUa = gO/lab + Oab + UJab , (6) 

where Va is the spatially totally projected covariant derivative operator orthogonal to u"-, Aa = Ua is the acceleration 
{Abv!' = 0), G is the expansion parameter, (Jab the shear {aab = (^(ab)j (^°'a = cFabu'' = 0) and ujab is the vorticity 
{(jJab = t^[a5]j i^abu'' = 0). Following the standard convention we will indicate the symmetrization over two indices of 
a tensor with round brackets and the antisymmetrization with square ones. 

In the u° frame, the Weyl or conformal curvature tensor Cabcd can be split into its Electric (Eab) and Magnetic 
{Hab) components respectively: 

Eab = Cacbd u'^ u'^ => E"' a = , Eab = E(^ab) > Eab = , (7) 
Hab = \ Vade C'^^bc -ff"a = , Hab = H(^ab) , Hab = . (8) 

In what follows we will use orthogonal projections of vectors and orthogonally projected symmetric trace-free part of 
tensors. They are defined as follows: 

Angle brackets may also be used to denote orthogonal projections of covariant time derivatives along u": 

v<'''> =h''bv'' , X^"''^ = [ - 5 /i"'' /led ] ^"'^ • (10) 

C. Energy-momentum tensors 

The choice of frame, i.e., choice of velocity field tt° and therefore the projection tensor hab allows one to obtain 
an irreducible decomposition of a generic Energy-Momentum Tensor (EMT), T^j*. The following unbarred quantities 
have been derived from the total EMT. quantities relating to the cfFcctivc fluids will be denoted with sub/super-scripts 
in order help to avoid confusion in later sections and to generalize to a multi-fluid system. 

Tab* = IJ-UaUb + phab + QaUb + QbUa + T^ab , (H) 

where fi is the total energy density and p is the total isotropic pressure of the fluid, qa represents the total energy 
flux, TTab is the total anisotropic pressure. Additionally, we have the following constraints 

QaU"' = , 7r"a = , TTab = T^{ab) , TTabU^ = ■ 

The various components of the total energy momentum tensor can be isolated in the following way: 



M = T*fu''u\ (12) 

P = l^fh'^'', (13) 

Qa = -Tlfu''h\, (14) 

ra6 = T'-' . (15) 



In a general fluid the pressure, energy density and entropy are related to each other by an equation of state 
p = p{ii, s) . A fluid is considered perfect if Qa and -Kab vanish, and barotropic if the entropy is a constant i.e. the 
equation of state reduces to p = p(/x). 
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D. Propagation and constraint equations 

Writing the Ricci and the Bianchi identities in terms of the 1 + 3 variables defined above, we obtain a set of evolution 



equations (here the 'curl' is defined as (JC)"'' = r?^''<° VcX%) : 

leg : rayQ - Vail" = - | 6^ + {uaii") - 2(7^ + 2a;2 _ i + 3p) , (16) 

leq : omrfotcj<"> - ^ r]"'"' VftUc = - § 6 + a^b a;^ (17) 

leq : sigdota^"^^ - VK''^ = - | 6 + u''^ - aK (J^^^ - J'^ - {E°-^ - \ tt"'') , (18) 

(£;<"''> + \ TT^-^^^) - {Hf'' + \ V<"g''> = - i (m + p) tr''^ - e (£:"^ + i tt'^'') (19) 

+ 3cr<"c(i^''^'=-g7r''>'=)-?i<'^9*> 
+ 77^''<" [ 2 tic if'^d + u;e (i^'^d + i Tr'-^d) ] , 

- rf''^'' [ 2 «e E''^ - \ a^'K - H% ] , 

leq : consl/i + Va?" = - 6 (/i + p) - 2 (?i„g") - ((T"67r''„) , (21) 

leq : cons2q^"-^ + V> + VfoTr"^ = - | 9 g° - a°6 g** - (/x + p) - W6 tt"^ - 7?"*"= Wfc gc , (22) 
and a set of constraints 

leq : divEVbiE'''' + i tt"^) - i V> + i 6 g» - i cr^b g*" - 3 if"^ - [ gm - ^uj^qc] = , (23) 

leq : divHVbH^^ + (m + p) + 3 (i;"*" - i tt"'') + r/«^^ [ i VftSe + {E\ + = , (24) 

leq : onuVw"^ - | V"e + t?"'"^ [ Vs^c + 2 Wc ] + = , (25) 

VaW" - {UaOj") = , (26) 

IhconstrH"^ + 2 w''^ + V<"w^> - (.a)"'' = , (27) 



that are completely equivalent to the Einstein equations. It is from these equations that we derive the general evolution 
equations for linear tensor perturbations. 

III. TENSOR PERTURBATION EQUATIONS 
A. The Background 

The equations presented in the previous section hold in any spacetime we may wish to analyze. However, in what 
follows wc will focus on the class of spacetimes that can be thought of as describing an " almost" Friedmann-Lemaitre- 
Robertson- Walker (FLRW) model, motivated by the fact that current observations suggest that the universe appears 
to deviate only slightly from homogeneity and isotropy. We can define a FLRW spacetime in terms of the variables 
above. Homogeneity and isotropy imply: 

a = a; = 0, Va/ = , (28) 
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where / is any scalar quantity; in particular 

VaAi = VaP = , Ua=0. (29) 

It follows that the governing equations for this background are 

e + ie2 + i(^ + 3p) = 0, (30) 

i? = 2[-ie2 + ^] , (31) 
A + e (m + p) = . (32) 

Now in order to describe small deviations from a FLRW spacetime we simply take all the quantities that are zero in 
the background as being first order, and retain in the equations (Eq. (??)-(??)) only the terms that are linear in these 
quantities, i.e. we drop all second order terms. This procedure corresponds to the linearization in the 1+3 covariant 
approach and it greatly simplifies the system of equations. In particular, the scalar, vector and tensor parts of the 
perturbations are decoupled, so that we are able to treat them separately. In what follows we will focus only on the 
tensor perturbations. 

B. The general linear tensor perturbation equations 

The 1+3 covariant description of gravitational waves in the context of cosmology has been considered by [2^. The 
linearized gravitational waves are described by the transverse and trace-free degrees of freedom once scalars have been 
switched off. Therefore, focusing only on tensor perturbations the necessary evolution equations are 

2 1 

O'afc + 2 O-afc + Eab - -TTab = , (33) 

Hat + Hate + iE)ab~^iTT)ab^O , (34) 

Eab + Eab Q - [H)ab + ^ (M + P) C^ab + 76 T^ab + 7;'^ab = , (35) 

2 5 2 

together with the conditions 

VbH^'' = 0, VbS"^ = 0, Hab = ia)ab- (36) 

Note that, since the linear tensor perturbations are frame- invariant, the structure of the equations does not depend 
on the choice of 4- velocity ,Ua. In the following, however, we shall choose the frame associated with standard matter 
{ua — u™). The motivation for such a choice is the fact that real observers are attached to galaxies and these galaxies 
follow the standard matter geodesies. Taking the time derivative of the above equations we obtain 

-5 /l 13\ 2 

<jab - ^^<J + 3 ^ + ( 9 ®^ + ^ 2^ ) = '^"fc + 3 ® '^ifc , (37) 

Hab - y^Hab + ^ e + ^ - 3p) Hab = (»afc + ^ © (. ^)ab , (38) 

Eab - ^^Eab + leEab+^ {Q^ ~ 3p) Eab ~ + p) {l + 3 C^) dab 



11-5 1 

■^T^ab - -^y^T^ab + g 6 TTafc + - (6^ - /i) TTab 



(39) 
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where = p/ (i and we have used the Raychaudhuri equation (Eq. (??)), the energy conservation equation (Eq. ((321) ') 
and the commutator identity 

{XU^{Xy,, + \{X)Q. (40) 

These equations generahze the tensor perturbation equations for an imperfect fluid that were derived in 31]. Once 
the form of the anisotropic pressure has been determined in Eq. ((37)) - ((39)) . the equations can be solved to give the 
evolution of tensor perturbations. As already noted in [2^ the presence of a term that contains the shear in Eq. p9|) 
makes this equation effectively third order, so that it is not possible to write down a closed wave equation for Eab- If 
TTah = 0, it is easy to show that for consistency, the solution for this field must also satisfy a wave equation because 
the shear is a solution of a wave equation and Eq. ([33]) holds. This will also be the case here because in our case 
TTah oc Gab and SO the anisotropic pressure will also satisfy a wave equation. 

Following standard harmonic analysis, Eq. ([37]) and Eq. (j38|) may be reduced to ordinary differential equations. It 
is standard ^24j to use trace-free symmetric tensor eigenfunctions of the spatial the Laplace-Beltrami operator defined 
by: 

V'Qafc = ^Qab , (41) 

where k = 27ra/A is the wavenumber and Qab — 0. Developing <Tab and Hab in terms of the Qab, Eq. (|37p and Eq. (j38p 
reduce to 

^^'^ + 1 e dC^) + {Iq' + 1^-Ip + 5) -f.^ = + 1 , (42) 

3 3 \ / 3 



3 2 



and Eq. ([33]) reads 



IV. GENERAL EQUATIONS FOR FOURTH ORDER GRAVITY 

The classical action for a fourth order theory of gravity is given by 

A = J d^x^[A + cnR + ciR^ + C2R^^R^"' +£m] , (45) 

where we have used the Gauss Bonnet theorem (33 | and Cm represents the matter contribution. In situations where 
the metric has a high degree of symmetry, this action can be further simplified. In particular, in the homogeneous 
and isotropic case the action for a general fourth order theory of gravity takes the form 

A = Jdx^ [/ (i?) + Cm] , (46) 

where Cm. represents the matter contribution. Such modifications to the linear Einstein-Hilbert action can typically 
arise in effective actions derived from higher-dimensional theories of gravity [20] . Varying the action with respect to 
the metric give the gravitational field equations: 

f'Rab~habf = {g\g\-gabg''')S,a + T2, (47) 
2 HV-g^m) 

matter and Sab = VaV6/'(-R). The trace of Eq. ^ gives: 

f'R-2f = + (48) 



where / = fiR), f — f'(R) = df(R)/dR, TlH — represents the stress energy tensor of standard 

'-g og^^ 
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where S = g Sab- The various components of Sab can be decomposed as 



Sab = /" 



VaVfoi? - VaRub - UaU^VciVbR) + RUaUb - R [V aUb - UaUb 



VaR^bR - R ( "^bRUa + Vai?Ub ) + R^UaUb 



S = f" (y^VcR + ii^VcR -R- or) + /"' (v^RVcR - R' 



(49) 



These equations reduce to the standard Einstein field equations when /(i?) = R. It is crucial for our purposes to be 
able to write Eq. (|T7|) in the form 



Gab 



rptot 

^ab 



rpR 

^ab^ 



(50) 



where T^^ = 



T, 



—r- and Tf? is defined as 
/' 



ab J, 



\{f- f'R) 9ab + {9\9\ - gabg'") Scd 



(51) 



The RHS of Eq. ([50)) represents two effective "fluids" : the curvature "fluid" (associated with T^) and the effective 
matter "fluid" (associated with T^). This step is important because it allows us to treat fourth order gravity as 
standard Einstein gravity in the presence of two "effective" fluids. This means that once the effective thermodynamics 
of these fluids has been studied, we can apply the covariant gauge invariant approach in the standard way. 

The conservation properties of these effective fluids are given by the Bianchi identities T^)J* When applied to the 
total stress energy tensor, these identities reveal that if standard matter is conserved, the total fluid is also conserved 
even though the curvature fluid may in general possess off-diagonal terms [HI, [H, [s^j ■ In other words, no matter how 



complicated the effective stress energy tensor T^^ is, it will always be divergence free if 
the single effective tensors, the Bianchi identities read 



rn:b 
b 



0. When applied to 



T. 



T, 



T. 



b 

R:b 
ab 



f" 

f f'^ ^ab 



ab 



JI2 '^ab 



TJl R'' 



(52) 
(53) 



with the last expression being a consequence of total energy-momentum conservation. It follows that the individual 
effective fluids are not conserved but exchange energy and momentum. 

It is worth noting here that even if the energy-momentum tensor associated with the effective matter source is not 
conserved, standard matter still follows the usual conservation equations T"l^ — 0. It is also important to stress that 
the fluids with and T^'^' defined above are effective and consequently can admit features that one would normally 
consider unphysical for a standard matter field. This means that all the thermodynamical quantities associated with 
the curvature defined below should be considered effective and not bounded by the usual constraints associated with 
matter fields. It is important to understand that this docs not compromise any of the thermodynamical properties of 
standard matter represented by the Lagrangian . 

In the matter frame the various components of the total energy momentum tensor, Eq. p2|) can be re- written 
in terms of the two effective fluids: 



TTob 



/' 



/' 

c 
/' 

^ab 
/' 



(54) 
(55) 
(56) 
(57) 
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where we assume that standard matter is a perfect fluid, i.e. = and tt^ = 0. The effective thermodynamical 
quantities for the curvature fluid are 



M = 77 



1 

/' 
1 

1 

1 

7' 



\{f Rf) + f'R + f'R^ + \<df"R - ^/"v^i? - ^/"'V'^i?v,i? - i/" i^.Vi? 

f'RVaR + /"Vai? - ^/"Vai? 



The twice contracted Bianchi Identities lead to evolution equations for /x™, , q^: 



J"R 



J"^aR 



(58) 
(59) 
(60) 
(61) 

(62) 
(63) 

(64) 



and a relation connecting the acceleration Ua to /i™ and follows from momentum conservation of standard matter: 

VV" = -(Ai™+p™)MV (65) 

Note that, as we have seen in the previous section the curvature fluid and the effective matter exchange energy and 
momentum. The decomposed interaction terms in Eq. and Eq. (|64p are given by /x™ ^ J,'i^ and ii^^-p^- 

It is easy to see that the curvature fluid is in general an imperfect fluid, i.e. has energy flux (g^) and anisotropic 
pressure {'Kab)- Since we are only interested in linear tensor perturbations, we need only be concerned with the tensor 
anisotropic pressure, which is proportional to the shear, Oab- We now present the second order evolution equations 
resulting from the standard harmonic analysis of Eq.s (l42 |) - (|44|) in the case of f(R) theories of gravity: 



1 /I 



fc2 
«2 



2 /' 



^Ir^U^'^^O, (66) 



-p 



n-2 



1 • f" 

-QR— 

3 /' 



--(/-/'i?)-i?2 



r 



- iiy \ = 0, 



(67) 



i?(^)--aW-0e + ii?^)aW. (68) 

For our purposes it will be particularly useful to consider these equations in the so-called Long Wavelength Limit. 
In this limit the wavenumber k is considered to be so small that the wavelength A — 2'Ka/k associated with it is 
much larger than the Hubble radius. Eq. (j4T|) then implies that all the Laplacians can be neglected and the spatial 
dependence of the perturbation variables can be factored out. 



V. TENSOR PERTURBATIONS IN R" GRAVITY 



To proceed, we must now fix our theory of gravity, i.e. we must choose the form of f{R)- We will consider a toy 
model (i?"-gravity) which is the simplest example of fourth order theory of gravity but exhibits many of the properties 
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of such theories. In this theory f{R) = X-R" a-nd the action reads 

A = J d^xV^[xR'' + CM] , 



(69) 



where % a the coupling constant with suitable dimensions and x = 1 for n = 1. li B. ^ 0, the field equations for this 
theory read 



— 1 ab I rpH 



nR 



where 



-1 ^ab 



^ab 



nR" 
in-l) 



R 
2n 



gab 



Ty-.cd T>-,c T>u 



{gcaOdb — gcdQab) 



(70) 

(71) 
(72) 



The FLRW dynamics of this model have been investigated in detail via a dynamical systems approach in 
where a complete phase space analysis was performed. This work demonstrated that for specific intervals of the 
parameter n there exist a set of initial conditions with non-zero measure for which the cosmic histories include a 
transient decelerated phase (during which large-scale structure can form) which evolves towards one with accelerated 
expansion. These transient almost Friedmann models existed for 0.28 < n < 1.35 in the case of a dust filled {w — 0) 
universe and for 0.31 ^ n < 1.29 in the case of a radiation filled {w — 1/3) universe. As we will discuss in later 
sections, these allowed intervals of n could be reduced significantly with future observations of gravitational waves. 
This model was also investigated as a possible explanation for the observed flatness of the rotation curves of spiral 
galaxies and the observed late times acceleration of the universe jl2| . The authors found a very good agreement 
between this model and observational data when n = 3.5. This is however at odds with the results of [ll[. Thus if one 
requires a transient decelerated phase (during which large-scale structure can form) and a solution to the dark matter 
and dark energy problem, the i?" model appears not to be viable. However the aim of this paper is to show that 
the study of tensor perturbations can in principle provide a strong constraint on the theory of gravity independent 
of existing cosmological data sets and consequently this work will provide a template for a more extensive study of 
tensor perturbations of f{R) cosmologies. 

In what follows we begin by analyzing the evolution of tensor perturbations in the absence of standard matter. We 
then consider the case of dust/radiation dominated evolution. Although we will give the full solutions, the discussion 
of the physics will be restricted to the long wavelength limit. 



A. The vacuum case 



We start by considering tensor perturbations in the absence of matter. This class of theories then admits the 
following exact solution 

ait)=a,t^, ^^(l^n)(2n-l)^ ^ ^ ^^3^ 



The expansion parameter is given by 



For the purposes of this paper we restrict our attention to expanding models. This requires q > 0, which in turn 
restricts the parameter n. In order to have an expanding background we require < ?i < 1/2 and 1 < n < 2 (we 
recover a static vacuum solution for n = 1/2, 1). We will only investigate models with values of n which satisfy 
the second inequality (since we wish to investigate models close to GR). The EOS of the total effective fluid in the 
background is then [j| 

p _ _1 {6if^7n-ll 
" M " 3(2n-l)(n-l)- ^^^^ 
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FIG. 1: The exponents of each mode of the solution Eq. [82] against n in the vacuum case. The black (grey) Une represents 
the growing (decaying) mode. The points represent the value of the exponents in the case of GR (n = 1). As n is increased 
the background expansion rate increase, resulting in a weaker growth rate for tensor perturbations. When the critical value of 
n 1.366 is reached no growing modes can be supported in this model. 



The EOS is singular and the pole's occur at n — 1/2 and n — I. Additionally, we have accelerated expansion, 
{w < —1/3) for n > (1 + V3)/2 « 1.366 and in the limit n — > oo we have w — > — 1. Substituting into Eq.'s (|66|) - (|68|) 
we obtain 



f(fc) 



3(1 - n)(4n- 3) 
{n-2)t ' 



.(fe) 



n{An-b){n- l)(8n-7) 
(n - 2)2t2 



a^^^ = 0, 



(76) 



(1 -n)(16n- 11) 
{n - 2)t 



2(6n2-8n + l)(n- l)(5n-4) 



H^^^ = 0, 



(77) 



{n - 2)t 

In the long wavelength/super-horizon limit (fc = 0), the above equations admit the following solutions 

- , „(4n-5) (S,.-7)(„~l) 

a^'') = Alt +A2t "-2 



H 



(fc) 



2(5n-4)(,i-l) 



,1.) ^ , ^ (n— l)(5n — 4) 8..^-16n + 9 

(n-2) 

However, the physical quantity of interest is the dimensionless expansion normalized shear E = a / H . 

,,\ ~ ~ (2ii-l)(4n-5) 

sc^) = Sit +S2r — ^^^^ — ^. 



(78) 

(79) 
(80) 
(81) 

(82) 



In Fig. [T] we have plotted the exponents of each mode of the solutions given above as a function of n in order to better 
see how the large scale behavior varies. The black (grey) line represents the growing (decaying) mode and the points 
represent the value of the exponents in the case of GR (n = 1). In the GR limit we recover a static vacuum model in 
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the background and always grows indicating that this model is unstable with respect to tensor perturbations. In 
the case of larger values of n, the Ei mode grows (decays) for n £ 1.366 {n ^ 1.366) and the S2 niode grows (decays) 
for n < 1.25 (n > 1.25). This is consistent with the background dynamics in that all perturbation modes are decaying 
when we have accelerated expansion {w < —1/3) in the background. 

For the sake of completeness we present the results of the general case (fc ^ 0). The solutions are given in terms of 
Bessel functions of the first and second kind (J and Y respectively). 



t '^(^^ 



(2n-l)(8n-9) 



AiJ (s,—j +A2Y is,— 
A3J [s,— ] +A4Y ( s, — 



12n^-22n + ll 

Alt 5f^J=5i 



rn2 + 2n-5 / kf 
J s, — 



12n''-22Ti + ll 

+A2t 



(2-n) 
+ 2n - 5 



kfJ 



kV 



(2-n) ^T' r 



kf 



+ kfY s + 1, 



kf 



where we have introduced the following parameters 

2n2 - 2n - 1 



n-2 



-1 + 



3(2n - 3) 
2(n- 2)r' 



The normalized shear S is now of the form 



, 3(4n^-6n+l) 

5](fe) = f — 2^^^^ — 



^ ke\ ~ ^/ kf 

EiJl s,— j + T,2Yls,— 



(83) 

(84) 

(85) 
(86) 

(87) 
(88) 



where both the Ej modes grow (decay) for n ^ 1.290 (n ^ 1.290). 



The fluid case 



We will now consider the case of tensor perturbations in the presence of matter which is described by a perfect 
fluid with barotropic EOS index, Wm- This class of theories then admits the following exact solution 



a(t) = 00^511+^, K = 0, 



The expansion parameter is given by 



e(t) = 



2n 



{1+Wm)t' 



(89) 



(90) 



As in the previous case we restrict our attention to expanding models. Additionally, we are mainly interested in 
the case where the perfect fluid describes dust {wm = 0) or radiation {wm = 1/3). This is due to the fact that these 
cases are the most relevant when considering GW detection via the CMB or direct detectors, e.g. LISA and BBO. To 
insure an expanding model we now require n > 0, provided Wm > — 1- 



1. The dust case 



We now investigate the evolution of tensor perturbations in the dust dominated era. The scale-factor is given by 

2n 



a(t) = aof, r 

' 1 1 1 ( EOS of the total effective fluid (dust and the effective curvature fluid) is then 

jn-l) 
w = —- -. 



(91) 



(92) 
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FIG. 2: The exponents of each mode of the solution for the normahzed shear against n in the dust and radiation dominated 
era's, (a) The left panel represents the exponents of the mode in the dust dominated era. (b) The right panel represents the 
exponents of the mode in the radiation dominated era. The black (grey) line represents the growing (decaying) mode. The 
points represent the value of the exponents in the case of GR (n = 1). 



The EOS is divergent for n = and we have accelerated expansion (w < —1/3) when n > 3/2. In the limit n —^ oo 
we have w — > — 1. Substituting into Eq.'s (|66p - (p5|) we obtain 



2(2n + 3) 
3t 



in -6) 



3^2 



= 0, 



(93) 



H 



(k) , 2(4n + 3) ^(,., 



3t 



2(2^2 + 5n- 3) 
3^2 



-It 



H^^^ = 0, 



(94) 



(n + 3) 

( 

it 



In the long wavelength limit (A: — 0), the above equations admit the following solutions 



(95) 



(96) 



H 



(k) _ 



-(r+2) 



2n) 



(97) 



The normalized shear is given by 



5n . 



(98) 



(99) 



The El mode is the decaying mode solution and is independent of the parameter n. This mode corresponds to the 
standard decaying mode found in GR. The E2 mode grows (decays) for n < 3/2 {n > 3/2) and reduces to the GR 
growing mode in the limit Ji 1. This is consistent with the background dynamics in that all perturbation modes 
are decaying when we have accelerated expansion {w < —1/3) in the background. 

In Fig. [2lja) we have plotted the exponents of each mode of the solutions given above as a function of n. The 
black (grey) lines represents the growing (decaying) mode and the points represent the value of the exponents in the 
case of GR (n — 1). For most of the values of n the perturbations grow slower in _R"-gravity than in GR. In fact 
only for n < 1, does the Ei mode grow with a rate faster than the usual t^/s. In the case of GR, there is always a 
growing tensor perturbation mode provided the background is not undergoing accelerated expansion. In the case of 
R" gravity, tensor perturbations grow at a slower rate, thus requiring a sufficiently decelerated expansion in order to 
support a growing mode. 
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Again, for the sake of completeness, in the general case (fc 7^ 0) the solutions are given in terms of Bessel functions 
of the first and second kind (J and Y respectively). 



^(fe) ^ ^-(2.+ l)/2 j^^j 



2r-3 /ct(i-'') 



2(r- 1)' (r- 1) 



2r-3 fct(i-^) 



2(r - 1)' (r - 1) 



(100) 



jj(k) ^ ^-(4,.+ l)/2 o J 



2r - 3 kt^^-'^ 



_2{r- 1)' (r- 1) 



BaY 



2r-3 H^i-'') 



2(r- 1)' (r- 1) 



(101) 



^(fe) ^ jg^^-(2.+3)/2| 2(2r 3) ^ 



2r - 3 /ct(i-'') 



+B2t 



-(2r+3)/2 



2(2r-3) 



2(r- 1)' (r- 1) 
2r - 3 kt'-^- 



1 



2(r- 1)' (r- 1) 



2(l-r)' (r-l)_ 
1 



2(1 -r)' (r- 1) 



(102) 



^. r/ie radiation case 

Next, we study the evolution of tensor perturbations in the radiation dominated era. The results of this section 
are especially relevant if one wishes to constrain f{R) models through their impact on the B-mode correlation on the 
CMB. The scale- factor goes as 



a(t) = aof ^ ^ 2 ■ 



The EOS of the total effective fluid (radiation and the effective curvature fluid) is then 

(3n - 4) 

w = . 

on 



(103) 



(104) 



The EOS is divergent for n — and we have accelerated expansion (w < —1/3) when n > 2. In the limit 71 —> 00 we 
have w ^ —1. Substituting into Eq.'s we obtain 



2t 



(4-n)(n-l) , ^2^_2r 



2^2 



a^'^) = 0, 



(105) 



rr{k) , 3n + 4 ■ (fe^ _ 

2t 



n{3n^2) 2^-2r 



t2 



+ kH- 



(106) 



t 

In the long wavelength limit (fc = 0), the above equations admit the following solutions 

E^^^ = Ci{2r-2)t-^'' + C2il-r)t'^'-^\ 



The normalized shear is given by 



(107) 

(108) 
(109) 
(110) 

(111) 



The El mode grows for < n < 2 and decays for n > 2. The S2 mode decays for the range < n < 2 and grows 
for n > 2. In Fig. [DJb) we have plotted the exponents of each mode of the solutions given above as a function of n. 
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The black (grey) lines represents the growing (decaying) mode and the points represent the value of the exponents 
in the case of GR (n = 1). For < n < 2 the Ei mode grows and the S2 mode decays. In the range n > 2 the 
modes change behavior in that the Ei mode decays and the E2 mode grows. Again, for most of the values of n the 
perturbations grow slower in i?"-gravity than in GR and only for < n < 1 and n > 4, does the modes grow with 
a rate faster than the usual linear growth. The most interesting feature of the solutions in the radiation dominated 
era is the possibility of growing modes even if the universe is in a state of accelerated expansion (n > 4). The impact 
of these modes on the CMB, could allow one to constrain deviations from GR. However, one should also analyze the 
evolution of perturbations on small scales. This analysis is beyond the scope of this paper and it is left to a future, 
more detailed investigation. 

In the general case {k ^ 0) the solutions are given in terms of Bessel functions of the first and second kind (J and 
Y respectively). 



3 fct(i 



2'(r-l) 



[r-l) 



(112) 



H(k) = i-(3r+l)/2 c J 



3 



2' (r-l) 



CiY 



3 



2' (r-l) 



(113) 



= Cit-(''+3)/2|2(r-l) J 



3 

2'(^ 

3 fci(i-'-) 
2'(^ 



- kt-^'Y 



5 yfci^i-'-) 
2'(^ 

5 



2' (r-l) 



(114) 



3. The generic large-scale case 

Finally, we study the evolution of large-scale [k — 0) tensor perturbations in the presence of a general barotropic 
fluid (thats is we will not fix Wm except to state that Wm > —1)- To insure an expanding model we now require n > 0, 
provided Wm > — 1. The scale-factor goes as 

2n 

a{t)=aoe r=—— -. (115) 

3 (l-f Wm) 

The EOS of the total effective fluid (radiation and the effective curvature fluid) is then 

{wm + l-n) 

w = . (116) 

n 

The EOS is divergent for n = and we have accelerated expansion {w < —1/3) when n > 3{w,n + l)/2. In the limit 
n ^ 00 we have w — > — 1. Substituting into Eq.'s ((66|) -([68 | we obtain 



„(fc) 10n-f6(l + w„0(l-n) ^^(fc) 2 (3 -f 3w„, - An) (nw^ - Wm - 1) ^(k) ^ g 
jj(k) I 14?^ + 6 (1 + w„0 (1 - n) ^(fc^ ^ 16n^ ~ Sn {wm + 1) + 6 {n - 1) {wm + 1) {w„i + 1 ~ 2n) ^(f^) ^ ^ 

E(k)- An + 3{l + Wm){l-n) ^^^)^ ^^^^^ 

+ Wm) t 



The solutions are then 



= T^itd-^-) + C2t('"-'-3''\ (120) 
= i:)i(n-2)t-2'-_^c'2(n-l-r)t(2«-3-3r)^ ^^22) 
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Region 1 
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-0.5 




w 
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FIG. 3: The range of parameters for which the modes grow or decay as a function of n and Wm The black thick line represents 
the change from growth to decay for the Ei mode. The grey thick line represents the change from growth to decay for the E2 
mode. The thin black line represents the GR case {n — 1). In region I, Ei decays and E2 grows. In region II, both modes 
grow. In region III, Ei grows and E2 decays. Finally in region IV, both modes decay. 

The normalized shear is given by 

^{k) ^ |,^^(2-2r) ^5.^^(2n-l-3r)_ ^^23) 

The El mode grows for n < 3{wm + l)/2 and decays for n > 3{wm + l)/2. The S2 mode decays for the range 
n < {wm + l)/2u;,„ and grows for n > (w,„ + l)/2wm- In Fig. [3] we have plotted the range of parameters for which 
the modes grow or decay as a function of n and Wm- This divides the parameter space into four regions. In region 
I, El decays and E2 grows. In region II, both modes grow. In region III, Ei grows and E2 decays. Finally in region 
IV, both modes decay. The most interesting feature of these solutions are those of region IV. As mentioned earlier 
this particular model was also investigated as a possible explanation for the observed flatness of the rotation curves 
of spiral galaxies and the observed late times acceleration of the universe [T^ . The authors found a good agreement 
between this model and observational data when n — 3.5 in the presence of dust (wm = 0). However, from our analysis 
we have found that such a choice of parameters ensures the absence of growing modes in the tensor perturbations. 
Therefore, if we wish to use this model as an explanation for Dark Matter, we can use gravitational wave detectors 
to severely constraint such theories. 

VI. CONCLUSIONS 

We have presented a mathematically well defined method of analyzing the evolution tensor perturbations of FLRW 
backgrounds in fourth order gravity, providing a general template for the study of linear gravitational waves in this 
context. The analysis is based on two important steps. Firstly, the recasting of the field equations for a generic fourth 
order theory of gravity into a form which is equivalent to GR, plus two effective fluids (the curvature "fluid" and 
the effective matter 'fluid"). Secondly, using 1+3 covariant approach, it is possible to derive the general equations 
describing the evolution of the cosmological perturbations of these models for a FLRW background. In this paper we 
have only dealt with the evolution of tensor perturbations, the evolution of scalar perturbations was been presented 
in ^iB'l and the vector perturbations will be presented elsewhere [ssj . Providing that one has a clear picture in mind 
of the effective nature of the ffuids involved, the approach above has the advantage of making the treatment of the 
perturbations physically clear and mathematically rigorous. 

Once the general perturbations equations were derived, we specialized them to the case of the i?"-gravity model. 
Using background solutions derived from an earlier dynamical systems analysis [ll| . we found exact solutions to the 
perturbation equations both in a vacuum and in the presence of matter (dust and radiation). We presented both the 
large-scale limit and full solutions, however, we restricted our discussions to the large-scale results. In Section IV Al 
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we studied the evolution of tensor perturbations in vacuum. The background solution proved to be unstable under 
tensor perturbations in the case of GR, where the background represents a static vacuum solution. In addition, for 
general values of n, the rate of growth of tensor perturbations is weaker than the GR case. 

In section IVBI we studied the evolution of tensor perturbations in the presence of matter. We first considered the 
case of the dust dominated era. For most choices of n (n > 1) the perturbations were found to grow at a slower rate 
in i?"-gravity than in GR and no growing mode could be supported for n > 3/2. 

Next, we studied the evolution of tensor perturbation in the radiation dominated era. Again, for most choices of n 
(n > 1) the perturbations were found to grow at a slower rate in i?"-gravity than in GR. However, it was found that 
their is always a growing mode present except for the special case of n = 2. This could have important consequences 
on the tensor perturbation spectrum, e.g. result in a tilt or running of the spectral index of the power spectrum. In 
this way the connection between the spectrum of tensor perturbations and the CMB polarization power spectrum 
offers an interesting independent way of testing for alternative gravity on cosmological scales. 

Finally, we studied the evolution of tensor perturbation in the presence of a generic fluid (wm > —1) in the large 
scale limit (fc = 0). We found that there is a range of the parameters Wm and n for which no growing modes are 
present (n > 3{wm + 1)/2 and n < {wm + l)/2wm). This corresponds the the choice of parameters as required to solve 
the Dark Matter problem in the work of |12l |. Thus providing an alternate method of constraining aforementioned 
theories of gravity. 

As in the case of the results found for the evolution of the scalar perturbations ^)3\ , the key question is how general 
these results are in terms of the form of the fourth order Lagrangian. Unfortunately this question is not easy to answer 
based only on the analysis presented above. The key point to consider would be the differences in the dynamics of 
the perturbations which, as we have seen, are very pronounced but more difficult to used because they depend largely 
on the features of the background. The important point, however, is that these differences do not necessarily imply a 
complete incompatibility with the data coming from the CMB and other observational constraints. Much more work 
will be needed before we can determine whether alternative gravity provides a viable alternative to standard General 
Relativity. 
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APPENDIX A: COVARIANT FORMALISM VERSUS BARDEEN'S FORMALISM 



As we have seen the covariant approach is a very useful framework for studying perturbations in alternative theories 
of gravity. However, since most work on cosmological perturbations is usually done using the Bardeen approach [2l| . 
we will give here a brief summary of how one can relate our quantities to the standard Bardeen quantities. A detailed 
analysis of the connection between these formalisms is given in . Here we limit ourselves to give the main results 
for tensor perturbations. 

In Bardeen's approach to perturbations of FLRW space-times, the metric gat is the fundamental object. If (jab is 
the background metric and gab — gab + 5 gab defines the metric perturbations 8gab in these coordinates. 
The perturbed metric can be written in the form 

ds^ = a2(7/){-(l + 2A)dry2 - 2B^dx''dr] + [(1 + 2ffi)7„^ + 2H^f^]dx''dx^} , (Al) 

where rj is the conformal time, and the spatial coordinates are left arbitrary. This spacetime can be foliated in 
3-hypersurfaces S characterized by constant conformal time rj and metric ^ab- 

The quantities A and Ba are respectively the perturbation in the lapse function (i.e. the ratio of the proper time 
distance and the coordinate time one between two constant time hypersurfaces) and in the shift vector (i.e. the rate 
of deviation of a constant space coordinate line from the normal line to a constant time hypersurface), Hj^ represents 
the amplitude of perturbation of a unit spatial volume and H^p is the amplitude of anisotropic distortion of each 
constant time hypersurface f25|. 

The minimal set of perturbation variables is completed by defining the fluctuations in the energy density: 



/i = /i + , 6 = , (A2) 

and the fluid velocity: 

M° = zZ" + Su" , 5u°' = u^v" , 5u^ = -u^A , (A3) 

together with the energy flux qa and the anisotropic pressure -Kab which are GI by themselves. 

These quantities are treated as 3-fields propagating on the background 3-geometry. With suitable choice of boundary 
conditions [36] , these quantities can be uniquely (but non-locally) decomposed into scalars, 3- vectors and 3-tensors: 



B^ = B\^+Bi, (A4) 

-S I ttTT 

T{a\f3) + ^Ta(i > 



where the slash indicates covariant differentiation with respect to the metric ^ab of S. In this way ^ abf — f\i3a~ 3^^/ 
and V^/ — /''''I'y is the Laplacian. The superscript 5 on a vector means it is solenoidal (sf " = 0), and TT tensors 

are transverse [H'^^^ \p = 0) and trace-free. 

On the base of IA4I and IA5[ it is standard to define scalar perturbations as those quantities which are 3-scalars, 
or are derived from a scalar through linear operations involving only the metric 7ab and its | derivative. Quantities 
derived from similar operations on solenoidal vectors and on TT tensors are dubbed vector and tensor perturbations. 
Scalar perturbations are relevant to matter clumping, i.e. correspond to density perturbations, while vector and tensor 
perturbations correspond to rotational perturbations and gravitational waves. 

Given the homogeneity and isotropy of the background, we can separate each variable into its time and spatial 
dependence using the method of harmonic decomposition. In the Bardeen approach the standard harmonic decompo- 
sition is performed using the eigenfunctions of the Laplace-Beltrami operator on 3-hypersurfaces of constant curvature 
E (i.e. on the homogeneous spatial sections of FLRW universes). In particular these harmonics are defined by 

y2y^fc) ^ _k2Y{k) ^ (A7) 

where F^^), rj^^ are the scalar, vector and tensor harmonics of order k. In this way one can decompose scalars, 
vectors and tensors as 

A = A{r])Y (A9) 

B„ = i3(0)(,7)ri")+i?«(ry)yji) , (AlO) 
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The key property of linear perturbation theory of FLRW space-times, arising from the unicity of the spHtting of (jA4p 
and (|A5p . is that in any vector and tensor equation the scalar, vector and tensor parts on each side are separately 
equal, i.e. the scalar, vector and tensor components of the equations decouple. 

All the quantities defined above can be decomposed in this way. However, before proceeding, one should note that 
the quantities A, Ba, , H^p,5,v" change their values under a change of correspondence between the perturbed 
"world" and the unperturbed background, i.e., under a gauge transformation. In order to have a gauge- invariant 
theory one has to look for combinations of these quantities which are gauge invariant. Bardeen constructed such GI 
variables to treat scalar and vector perturbations (2l| . The quantities which are relevant to our analysis, tTq,/? and 

^Tap (^^^ harmonically decomposed object Hip) are already GI. 

The variables covariantly defined in the main text are, by themselves, exact quantities (defined in any space-time) 
and are GI by themselves, therefore, to first order, we can express them as linear combinations of Bardeen's GI 
variables. In 2j| these expansions are given in full generality. Here we will limit ourselves to a few examples, giving 
only the tensor contributions and refer the reader to [131 for details. 

The tensor part of the shear, trace-free part of the 3-Ricci tensor, the electric and magnetic parts of the Weyl tensor 
are given by 



(2)/T^(2) 



af3 



(A12) 



E, 



1 r 
2 



H. 



(2)// 



2K) H. 



(2) 



(2) 



(A13) 
(AM) 



H. 



a/3 



„-2 rr(2)/T^(2)T.|5 



(A15) 



where the prime denotes derivative with respect to the conformal time 77. The relations above can be used to give 
an intrinsic physical and geometrical meaning to Bardeen's variables, and also to recover his equations. For example, 
combining our linearized expression for the trace-free part of the 3-Ricci tensor 



(3)7^ 



e 



Q/3 



E, 



1 



Q/3 



(A16) 



with the above expressions (Eq. IA12IIA14)) gives Bardeen's expression for the transverse and trace-free metric pertur- 
bation evolution equation 



(2)// 



2K) H, 



(2) 



(A17) 



where tt is the harmonically decomposed anisotropic pressure. Substituting for tt using Eg. 1611 and Eq. IA121 we find 
the general evolution equation for tensor perturbations in fourth order gravity theories to be 



(2)// 



d'f fdf 



2^- 

a ' 9i?2 \dR 



^ R' 



H, 



(2)/ 



2K) = 0, 



(A18) 



where primes denote differentiation with respect to conformal time throughout this appendix. 
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